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General relativity predicts the gravitational radiation signatures of mergers of compact binaries, 
such as coalescing binary black hole systems. Derivations of waveform predictions for such systems 
are required for optimal scientific analysis of observational gravitational wave data, and have so far 
been achieved primarily with the aid of the post-Newtonian (PN) approximation. The quality of this 
treatment is unclear, however, for the important late inspiral portion. We derive late-inspiral wave- 
forms via a complementary approach, direct numerical simulation of Einstein’s equations, which has 
recently matured sufficiently for such applications. We compare waveform phasing from simulations 
covering the last ~ 14 cycles of gravitational radiation from an equal-mass binary system of non- 
spinning black holes with corresponding 3PN and 3.5PN waveforms. We find phasing agreement 
consistent with internal error estimates based in either approach, at the level of one radian over 
~ 10 cycles. The result suggests that PN waveforms for this system are effective roughly until the 
system reaches its last stable orbit just prior to the final merger. 

PACS numbers: 04.25.Dm, 04.25.Nx 04.30.-w 04.30.Db, 95.30.Sf, 97.60.Lf 


Compact astrophysical binaries spiral together due to 
the emission of gravitational radiation. Calculating the 
dynamics of these systems, and the corresponding gravi- 
tational waveforms, has been a central problem in general 
relativity for several decades. With first-generation inter- 
ferometric detectors such as LIGO,. VIRGO, and GE0600 
now operating, and development moving forward on the 
space-based LISA mission, the need for accurate and re- 
liable waveforms for gravitational wave data analysis has 
become urgent. 

Post-Newtonian (PN) methods, based on expansions 
in the parameter e ~ v/c, have been the major analytic 
tool used to calculate the system dynamics and wave- 
forms during the early part of the inspiral, when the 
binary components are relatively widely separated and 
thus have a small orbital frequency [1]. Currently, grav- 
itational wave data analysis for binary inspiral relies on 
waveforms derived from PN methods [2]. The current 
predicted orbital phase is available up to 0(e 7 ), which 
is referred to as 3.5PN order. However the convergence 
properties of the PN sequence are not well understood, 
and it is not yet clear how well PN predictions work late 
in the inspiral when frequencies are high. 

Numerical relativity, in which the full set of Einstein’s 
equations is solved on a computer, is needed to handle 
the final stages of the binary evolution, when the com- 
ponents inspiral rapidly and merge. Recently, there has 
been dramatic progress in the use of numerical relativity 
to simulate the final inspiral and merger of black holes 
[3-11]. These breakthroughs have allowed numerical sim- 
ulations with increasingly wider initial separations, pro- 
ducing longer wavetrains. Linking such simulations with 
the PN calculations and comparing their waveforms in 
the late inspiral regime is a pressing concern of gravita- 


tional wave data analysis. 

We have carried out numerical simulations of a merg- 
ing equal-mass, nonspinning black-hole binary of suffi- 
cient duration to conduct meaningful comparisons with 
the PN approximation. The black holes start on nearly 
circular orbiting trajectories ~ 1200M before merger, 
where M refers to the mass that the system would have 
had when the black holes were still far apart, before ra- 
diative losses were significant. M is related to time by 
M = 5x10 ~ 6 s(M/Mq). In this Letter, we quantitatively 
compare crucial phasing information in the our numeri- 
cal simulation waveforms with phasing in PN waveforms, 
finding striking agreement. 

The numerical' simulations were performed using the 
moving puncture method [5, 6, 12]. We use fourth-order 
Runge-Kutta time integration, fourth-order accurate fi- 
nite spatial differencing, and second-order-accurate ini- 
tial data. Adaptive mesh refinement is used to resolve 
both the dynamics near the black holes and the propaga- 
tion of the gravitational waves [8]. We performed physi- 
cally equivalent runs at three different maximum resolu- 
tions: low (3M/64), medium (3M/80), and high (M/32). 
We find fourth-order convergence of the Hamiltonian con- 
straint, and better than second-order convergence of the 
momentum constraints during the runs. 

The simulations begin at an angular gravitational wave 
frequency w ~ 0.051 M _1 . The frequency then sweeps 
upward through roughly an order of magnitude while 
the black holes undergo ~ 7 orbits, thus producing ~ 14 
gravitational wave cycles before merger. For such long- 
lasting simulations, the primary consideration in provid- 
ing a realistic initial data model is to set up the orbiting 
black holes with minimal eccentricity, as gravitating bi- 
nary systems of comparable- mass objects are expected to 
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FIG. 1: Gravitational strain waveforms from the merger of 
equal-mass Schwarzschild black holes. The solid curve is the 
waveform from the high resolution numerical simulation, and 
the dashed curve is a PN waveform with 3.5PN order phas- 
ing and 2.5PN order amplitude accuracy. Time t = 0 is the 
moment of peak radiation amplitude in the simulation. 


circularize rapidly through the emission of gravitational 
radiation. We have selected an initial black hole config- 
uration with relatively little eccentricity of less than one 
percent, as measured below. 

Fig. 1 shows a comparison of gravitational wave strain 
generated by our highest-resolution numerical run and 
that predicted by PN with 3.5PN phasing and 2.5PN am- 
plitude accuracy. The waves are based on the dominant 
l = 2, to = 2 spin- weighted spherical harmonic of the 
radiation, and represent an observation made on the sys- 
tem’s equatorial plane, where only one polarization com- 
ponent contributes to the measured strain. The initial 
phase and initial time of the waves have been adjusted 
so that the frequency and phase for each waveform agree 
early in the simulation, at t = — 1000M. We will quanti- 
tatively study the evident phase agreement below. 

To conduct comparisons with PN calculations, we need 
to extract an instantaneous gauge-invariant polarization 
phase and angular frequency w from our simulations. 
These are derived from the first time-derivative of the 
gravitational wave strain, which is a robust quantity in 
the numerical data. This frequency corresponds to the 
sweep rate of the polarization angle of the circularly po- 
larized gravitational wave that can be observed on the 
system’s rotation axis. 

We define eccentricity as deviation from an underlying 
smooth, secular trend. We obtain a monotonic “secular” 
frequency-time relation by modeling the frequency u> as 
a fourth-order monotonic polynomial w c (f), plus an ec- 
centric modulation in the waveform angular frequency, 
du>(t), of the form du>(t) = oj(t) — w c (t) = Asin($(t)), 
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FIG. 2: Pit to estimate eccentricity in the high resolution 
run. The dashed (green) curve is the numerically observed 
frequency. The solid (red) curve is a fit to a secular increase 
with time with a sinusoidal perturbation caused by the ec- 
centricity. The dotted (blue) curve is the “de-eccentrized” 
secular trend <n c . 


where <h(f) is a quadratic function of time. The fitting 
procedure is illustrated in Fig. 2, which shows the orig- 
inal numerical frequency w(t), the fit, and tu c (f). Fit- 
ting this data yields A — 8(±1) x 10~ 4 . For Keple- 
rian systems, conserved angular momentum is propor- 
tional to r 2 u> so the (radial) eccentricity corresponds 
to half the fractional amplitude of frequency modula- 
tions: e — A/(2uj). In our case the eccentricity starts 
near 0.008, decreasing by a factor of three by the time 
uj c M ~ 0.15. We will compare our simulation with 
non-eccentric PN calculations, with the expectation that 
small eccentricities have minimal effect on the important 
underlying secular trend in the rate at which frequency 
sweeps up approaching merger. 

The phasing of the waveform is critical for gravita- 
tional wave observation. For data analysis, the optimal 
methods for both detection and parameter estimation 
rely on matched filtering, which employs a weighted in- 
ner product that can be expressed in Fourier space as 
< h, s >= / (h* (f)s(f) + h(f)s*(f))/Sn(f)df, where h is 
the template being used, s is the signal being analyzed, 
and S n is the one-sided power spectral density of the de- 
tector’s noise [13]. A template that maximizes < h,s > 
will provide an optimal filter. While a template with 
time-varying amplitude can emphasize some portions of 
a signal and not others, the crucial factor is the relative 
phasing of the template and signal. The inner product 
will cease to accumulate in sweeping through frequency 
if the template and the signal evolve to be out of phase 
with each other by more than a half-cycle, decreasing the 
effectiveness of the procedure. 

Our key objective is to compare phasing between nu- 
merical and PN waveforms. To avoid issues with timing 
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alignment, we will compare phases as a function of po- 
larization frequency, which corresponds to twice the or- 
bital frequency in the PN case. For circular inspiral this 
frequency should grow monotonically in time, with the 
frequency oj c providing a physical reference of the “hard- 
ness” of the tightening binary. 

Circular inspiral phasing information is typically de- 
rived in PN theory by imposing an energy balance re- 
lation to deduce the rate at which w c evolves from the 
radiation rate at a specified value of w c [1]. Though not 
strictly derived in the PN context, this physically sensi- 
ble condition currently allows the determination of the 
chirp rate w c (tn e ) (or something equivalent) up to 3.5PN 
order [14]. From such a relation information about phase 
and time are determined by integrating d4>/du> c = w e /w c 
and dt/div c = l/w c . The phasing information can be rep- 
resented by any one of several relations between phase, 
frequency and time. Various approaches take the PN- 
expanded representation of one of these relations as the 
PN “result” for waveform phasing [1, 3, 14]. We consider 
both the PN expansion of the chirp rate, numerically in- 
tegrated as needed, as well as the PN expansion of the 
phase. 

For the purpose of comparison with our numerical sim- 
ulations, we invert the monotonic function u c (t) to obtain 
the phase as a function of frequency: = <p(t(ui c ))- 

Note that the effect of eccentricity is not removed from 
cj>, though the “circularized” frequency co c does provide 
the abscissa according to which phases are compared in 
the different treatments. , 

Fig. 3 shows the wave phases; here the phases are ad- 
justed by addition of a constant so that each phase van- 
ishes at u c M = 0.054, corresponding to the time 1000M 
before the radiation merger peak in the high resolution 
numerical simulation. The sequence of numerical results 
converges at fourth-order (verified in Fig. 4), allowing us 
to project by Richardson extrapolation to a fifth-order 
accurate result (thick solid line in Fig.3). 

In Fig. 3 we compare the numerical simulation results 
with two versions of PN-phasing, based on either the PN 
expansion of the phase or on the numerically integrated 
PN expansion of the chirp rate. We show each of these 
at 3PN and 3.5PN order. While the 3PN phase seems 
to agree closely with the extrapolated numerical phase, 
the 3.5PN phase moves away, and develops an unphysi- 
cal negative slope after u c M ~ 0T5. The agreement of 
the extrapolated numerical result with the integrated PN 
chirp rate improves with PN order, with the 3.5PN result 
showing striking agreement up to about w c M ~ 0.15. 

We look more quantitatively at the differences among 
the phases in Fig. 4. The solid curve shows the differ- 
ence between our medium and high resolution results, 
while the dotted curve shows the difference between our 
low and medium resolution results, scaled such that for 
fourth-order convergence the curves should superpose. 
This is indeed demonstrated to good approximation. A 



FIG. 3: Gravitational wave phase, in radians, for numeri- 
cal and two versions of PN waveforms. The sequence of nu- 
merical simulation results approaches the solid curve, which 
agrees well with the phase obtained by numerically integrat- 
ing the 3.5PN expansion of the chirp rate w c (u> c ). It also 
agrees with the PN expansion of 4>{ls c ) at 3PN order though 
the 3.5PN order curve reaches an unphysical maximum value 
near ui c M ~ 0.2. 



FIG. 4: Gravitational wave phase error estimates. Differ- 
ences between phasing from the integrated 3.5PN chirp rate 
and Richardson extrapolation from the numerical simulations 
(solid curve) are small, and are consistent with internal error 
estimates for the numerical simulation results and the PN 
sequence. 

good estimate for the error of the phase in the high resolu- 
tion run is given by its difference from the phase obtained 
by Richardson extrapolation; this comes out to ~ 93% 
of the med-high (solid) curve shown in the figure. Note 
that the cumulative errors in the numerically-generated 
waveforms accrue primarily before the orbital frequency 
w c M = 0.1. (This makes sense generally since the simu- 
lations spend much more time at lower frequencies.) 
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Without monotonic convergence between the 2PN, 
2.5PN and 3PN at the frequencies considered here, it is 
difficult to estimate errors in the PN phase. Nonetheless 
we tentatively take the difference between the integrated 
3PN and 3.5PN chirp rates, shown by the dashed line in 
Fig. 4, as a measure of PN errors. 

The trend in the slope of these error curves indicates 
the rate at which phase error accumulates, as indepen- 
dently estimated within each approach. Fig. 4 sug- 
gests that phasing errors for our high resolution simu- 
lation accumulate more quickly than PN phasing errors 
for ui c M < 0.08 ( t < —300 M), with the numerical sim- 
ulation phasing being more accurate than PN at higher 
frequencies. In both cases, the phase error accumulates 
to roughly two radians by oj c M ~ 0.2 as the black holes 
begin to plunge together. 

We now address the central objective of this letter, 
a quantitative comparison of numerical and PN phasing 
results. We compare the Richardson-extrapolated phase, 
our best simulation estimate, with the integrated 3.5PN 
chirp rate result. The difference is shown in the solid 
curve of Fig. 4. Note that the phase differences be- 
tween the PN and numerical results (solid curve) accu- 
mulate less rapidly than the estimated errors in either 
the high resolution simulation or the 3.5PN results over 
much of the frequency range shown. This indicates that 
the numerical and post-Newtonian predictions appear to 
be converging to a common answer for waveform phase, 
in a frequency regime where the validity of the PN pre- 
dictions could not previously be assessed. The phase dif- 
ference among the best estimates of the numerical and 
PN approaches begins to grow more significantly after 
about uj c M ~ 0.15 ( t ~ —70 M) occurring about a cycle 
before the estimated end of the binary’s last stable orbit. 

Through the frequency range 0.054 < u> c M < 0.15 
the net phase difference, measured against frequency, 
amounts to less than one radian, a level of error which 
would be tolerable in many gravitational wave data anal- 
ysis applications. For most of the range studied here, the 
frequency-based phase differences between our high reso- 
lution simulation and the integrated 3.5PN chirp rate re- 
sult seem larger than time-based phase differences, which 
are evident in Fig. 1. This is because the sweep rate dif- 
ferences have not yet grown to a significant phase differ- 
ence; the difference would become evident if the wave- 
form was continued farther. At the higher frequency end 
of the waveform, the time-based differences would appear 
larger, unless the waveform was aligned to agree at the 
end of the waveform. This sensitivity to time-alignment 
makes frequency-based phase comparisons a more reli- 
able indicator of phasing differences. 

Our results provide a crucial cross-validation of PN 
waveforms from the late inspiral of binary merger, with 
results of new long-lasting numerical simulations. These 


simulations have sufficient accuracy to provide a mean- 
ingful comparison with PN waveforms over the last t ~ 
1000M of the coalescence, specifically addressing a bi- 
nary system of equal-mass non-spinning black holes. We 
find phase agreement consistent with internal phase-error 
estimates conducted in each approach, indicating that 
phase accuracies within a few radians are now achievable 
for this part of the coalescence waveform. 

We emphasize, however, that there is still much impor- 
tant work to be done in improving and futher assessing 
PN and numerical simulation waveforms. Certainly we 
have only addressed one case in a large parameter space 
of potential binaries, which will inevitably include sys- 
tems, such as rapidly precessing unequal-mass spinning 
binaries, which are harder to treat with present PN and 
numerical techniques. With either approach, even for our 
simple case, a non-negligible amount of phasing error ac- 
cumulates over the range studied, and more will have ac- 
cumulated at lower-frequencies addressible through the 
PN approximation. We expect continuing developments 
in numerical simulations and the pursuit of higher-order 
PN treatments to be crucial for developing a refined un- 
derstanding of coalescence waveforms, which will be cru- 
cial in some data analysis applications for gravitational 
wave observations. 
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